Drifting solitary waves in a reaction-diffusion medium with differential advection 
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Propagation of solitary waves in the presence of autocatalysis, diffusion, and symmetry breaking (differential) 
advection, is being studied. The focus is on drifting (propagating with advection) pulses that form via a convec- 
tive instability at lower reaction rates of the autocatalytic activator, i.e. the advective flow overcomes the fast 
excitation and induces a drifting fluid type behavior. Using spatial dynamics analysis of a minimal case model, 
we present the properties and the organization of such pulses. The insights underly a general understanding of 
localized transport in simple reaction-diffusion-advection models and thus provide a background to potential 
chemical and biological applications. 
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Solitary waves are prominent generic solutions to reaction- 
diffusion (RD) systems and basic to many applied science dis- 
ciplines ill] . In one spatial dimension, these spatially localized 
propagating pulses are qualitatively described by a fast exci- 
tation (leading front) from a rest state followed by a slow re- 
covery (rear front) to the same uniform state [1]. Thus, in 
isotropic RD media a single symmetric supra-threshold lo- 
calized perturbation results in counter-propagating pulses or 
wave trains 10]. 

However, in chemical and biological media transport can be 
facilitated by both diffusion and advection and thus excitation 
properties of solitary waves can be subjected to convective 
instabilities Nevertheless, theoretical foundations of soli- 
tary waves in the presence of a symmetry breaking advective 
transport, i.e. dynamics in a differential reaction-diffusion- 
advection (RDA) media, are yet to be established. Among a 
few reported examples, it was only shown both experimen- 
tally and numerically (with no underlying theoretical basis), 
that excitable pulses can persist in RDA with a propagation 
direction against the advective field [4], i.e. upstream. 

In this Letter we analyze an RDA case model and demon- 
strate that under low reaction rates, solitary waves may be- 
come convectively unstable and thus drift (see Fig.[T]), i.e. the 
slow recovery becomes a leading front. We reveal the regions 
and the properties of such drifting pulses and show that the 
phenomenon underlies a competition between a local kinetics 
of the activator and a differential advection. Our methods in- 
clude a bifurcation theory of coexisting spatial solutions (lin- 
ear analysis and numerical continuations) coupled to temporal 
stability; all the results well agree with direct numerical inte- 
grations. At the end, we discuss the potential applicability of 
our findings to chemical and biological media. 

Model setup - We start with a minimal RDA model that 
incorporates local kinetics of activator v{x^t) and inhibitor 
u{x^t) type: 

ut + sux = f{u,v) - u, (1) 
Levt^svx = Bf(u,v) - av ^ Pe~^Vxx. 

These dimensionless equations describe a membrane (or 
cross-flow) reactor, with continuous feeding and cooling in 



which an exothermic reaction takes place f{u^v) = Da{l — 
u) exp [7^/(7 + v)] [6]. Eq. ([T]) admits a uniform rest state 
{u^v) = {uq^vq = Buo/a), where uq obtained via Da = 
1^0 (1 — '^o)""'^ exp [—'yuo/{ja/B -\- uq)]. In what follows, we 
set Pe = 15, s = 1, a = 4, 7 = 10000, and use Le, Da and 
B as control parameters allowed to vary; parameter definitions 
are given in [7] . 

A standard linear stability analysis to periodic perturba- 
tions, shows that the uniform states {uq^vq) may loose sta- 
bility to two finite wavenumber Hopf instabilities, Da"^, that 
emerge from (P^, Da^), as shown in Fig.[Tl the instabili- 
ties are of a drifting type, i.e. in direction of advection. This 
anomaly arises due to the broken reflection symmetry of left- 
right traveling waves that is preserved in RD systems. We 
note that traveling waves TW~ bifurcate (nonlinearly) sub- 
critically from Da~ while traveling waves TW^ bifurcate su- 
percritically from Da^ [5]. While the region Da~ < Da < 
Da^ is linearly unstable, under certain conditions stationary 
periodic (SP) solutions may also develop. The criterion for 
SP states is zero of the real and the imaginary parts in the 
dispersion relation (for a finite wavenumber), identifying zero 
speed [5] (see dotted line in Fig.[T]). 

Here, our interest is in the affect of a differential advec- 
tion (Le) and the local kinetics (P, Da) on the organization 
of drifting solitary waves. We also consider large domains in 
which pulse behavior is not affected by the type of boundary 
conditions (periodic, no-flux or mixed) and also not interested 
in the regimes in which nonuniform steady state patterns may 
form, for details on the affect of boundary conditions see (3] . 

Propagation of solitary waves.- To reveal the propagation 
properties and the regimes of solitary waves (see Fig. [U, we 
look at the steady state version of ([B in a comoving frame, 
^ = x-ct [5]: 

= {s — c)~^ [Daf{u, v) — u] ^ = (2) 
— Pe [{s — cLe) w — BDafiu^ v) + av] . 

The advantage is that existence of nonuniform states can be 
now analyzed via spatial dynamics methods, i.e. where space 
is viewed as a time-like variable. Thus, solitary waves [in the 
context of ([T])] become in Q asymmetric homoclinic orbits 
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FIG. 1: (color online) Top panel: Regions of excitable (thick dashed 
line) and drifting (thin dashed line) solitary waves (pulses) in a reac- 
tion rate parameter space (B, Da) at Le — 100; the thin dashed line 
implies zero velocity of a pulse. The solid line mark the onsets of fi- 
nite wavenumber instabilities of travelins waves, TW^. The dotted 
line marks the criterion for stationary periodic (5'P) solutions (see 
text for details). The (•), marks the leftmost limits of homoclinic 
orbits (B*, Da*) ^ (8.76, 0.155) and asymmetric finite wavenum- 
ber Hopf bifurcation {B^ ,Da^) ^ (10,0.36) while (■) marks 
the leftmost limit of excitable pulses {Bo,Da'') ^ (10.35, 0.155). 
Bottom panel: Space-time plots at S = 10.6 and (a) Da = 0.24, (b) 
Da = 0.15, and (c) Da = 0.08. The plots show v{x,t) resulting 
from integration of Eq. ([T]) with no-flux boundary conditions, where 
X G [0, 10] and t G [0, 2265]; we used a top-hat initial condition 
embedded in {uo,vo) at the respective Da values. 

(HO) and TW^ (which will be also discussed) correspond 
to periodic orbits undergoing Hopf bifurcations at Da^ (with 
a proper c). In the following all these solutions will be com- 
puted numerically using a continuation package AUTO Isl], 
where the speed c is obtained as a nonlinear eigenvalue prob- 
lem. Then temporal stability of such steady states will be cal- 
culated employing an eigenvalue problem via Eq. ([T]) in a co- 
moving frame. 

In Eig. [2top panel), we present the branches of HO at 
Da = Da"" 0.155 (a horizontal cut in top panel in Eig.[T]), 
resulting via a simultaneous variation of (5, c). B = B* 
8.76 identifies a fold, where the stable branch corresponds to 
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FIG. 2: (color online) Bifurcation diagrams showing the branches 
of homoclinic solutions as a function of B at Da* 0.155 and 
Le — 100 (top panel) and Le — 1 (bottom panel). The branches 
are plotted in terms of the propagation speed and the maximal value 
of ^;(^) (top inset); solid lines indicate linear stability. Bottom insets 
show profiles ofv{^) at the two folds, marked by (•). The branches 
were obtained via integration of Eq. ^ while the stable portions 
of each branch coincide with solutions obtained by integration of 
Eq. IT]); here the periodic domain is L = 24 but results are identical 
also on larger domains. 



large amplitude HO (see inset). The drifting pulses exist for 
B* < B < Bq 10.35 since both branches have positive 
speeds, and have similar profiles (see inset) as the standard 
excitable pulses. Namely, drifting pulses propagate in the di- 
rection of the advection (downstream, c > 0) where the lead- 
ing front is the oscillatory tail that was a trailing tail above 
B = Bq, for excitable pulses (upstream, c < 0). 

Drifting pulses are expected at low reaction rate regimes of 
the activator, represented in ([T]) by dimensionless rate constant 
(Da) and exothermocity (B). Under such conditions the exci- 
tation of nearest neighbors is suppressed due to the advective 
flow (a nonlinear convective instability) and thus the pulse af- 
ter speed reversal is no longer excitable since the leading front 
now develops from the rest state as a small amplitude pertur- 
bation. This scenario changes once the differential advection 
is eliminated (Le = 1), in this case a typical RD behavior is 
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FIG. 3: (color online) Top panel: Schematic representation of typical 
eigenvalue configurations about the uniform state (i^o, '^^o, 0) corre- 
sponding to a saddle if B < Bb and a saddle-focus if B > Bh, where 
5b ~ 9.1 is the Belyakov point. Bottom panel: Typical dispersion 
relations that associated with the respective eigenvalues computed 
starting from the stable homoclinic orbits (see top panel). Parame- 
ters as in the top panel of Fig [2] 
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FIG. 4: (color online) Bifurcation diagram showing the branches of 
traveling waves {TW~) as a function of B in terms of speed and 
the maximal value of v{^) (in the inset), B — 10.4 (dark line) 
and B = 10.2, where Da~ 0.29, /cc ^ 3.2, c ^ 0.0054 and 
Da~ 0.31, /cc ^ 3.355, c 0.0053, respectively. Solid lines 
imply linear stability to long wave lengths perturbations |9], while 
(•) marks the respective onsets of the linear finite wavenumber Hopf 
bifurcation to TW~ . Integration details as in the top panel of Fig. [51 
but on distinct periodic domains. 



restored. While the c = line for Le = 1, in the (5, Da) 
plane doesn't change, we show in Fig. ^bottom panel) that 
near the fold only a negative velocity region forms, i.e. a stan- 
dard excitable pulses are being restored (stability of the pulses 
does not play a qualitative role). 

Nevertheless, drifting pulses in presence of a differential 
advection inherit the properties of excitable pulses, as demon- 
strated by monotonic and nonmonotonic dispersion relations 
in Fig. [3l The latter are important characteristics of organi- 
zation and interaction of solitary waves |10] and are distin- 
guished here around B = B^ 9.1, a so called Belyakov 
point I 111 . At this point and with an appropriate speed, the 
spatial eigenvalues [of Eq. Q] correspond to one positive real 
(associated with ^ — oo) and a degenerate pair of nega- 
tive reals (associated with ^ oo). Below B^, the degen- 
eracy is removed but the eigenvalues remain negative reals 
(a saddle) while above B^ they become complex conjugated 
corresponding to a saddle-focus (a Shil'nikov type HO 1 12]), 
marked by ( x ) in top panel of Fig. [3l Importantly, such an in- 
terchange of eigenvalues implies a transition from monotonic 
to oscillatory dispersion relation [Fig.^bottom panel)] and a 
monotonic (in space) approach of the HO to the fixed point as 
^ — > ±oo, which implies coexistence of bounded-pulse states 
for 5 > Bi, LIOJ. 

Organization of drifting states.- A standard theory of soli- 
tary waves qualitatively predicts an organization of HO to be 
accompanied by periodic solutions [12]. Here the dispersion 
relations obtained at B < B^ [Fig.^bottom panel)], indeed 
imply existence of periodic orbits although the uniform state is 
linearly stable. These periodic solutions are in fact TW~ that 



bifurcate subcritically from the locus of points Da = Da~ 
for 5 > B^ [with distinct critical wavenumbers and speeds 
obtained from the linear analysis of Eq. ([T])], as shown by two 
examples in Fig.|4l Notably, there are infinite number of such 
TW- families. Unlike the HO, stability of TW' solutions 
do depend on domain size [9] . 

The organization of all drifting nonuniform solutions can 
be understood by varying Da at two representative B val- 
ues. Fig. Oa), shows a bifurcation diagram of nonuniform 
solutions at B 10.4: while TW"^ propagate downstream. 
The single pulse HO branch ends at the two rightmost ends 
(marked by dots), at which the profiles take the form of ho- 
moclinic tails (see bottom inset) [13]. Due to the proximity 
to the subcritical onset of TW~ at Da~ , the two rightmost 
ends ever approach each other as domain (L) is increased, and 
consequently, they inherit the propagation direction of the top 
and the bottom branches of TW~ as discussed in |5]. As 
B is decreased below B^ the HO and the TW~ solutions 
organize in isolas and parts of their stability regions overlap 
[Fig. Ob)], implying sensitivity to initial perturbations. Note 
that the oscillations of the right tail in the profile had decreased 
(see bottom inset), which is consistent with the approach to- 
wards Belyakov point (B = B})). 

Conclusions and prospects- We have showed that solitary 
waves can propagate bidirectionally (without changing their 
shape) due to a competition between activator autocatalysis 
and a symmetry breaking advection. Consequently, we distin- 
guish between excitable (upstream or against advection) and 
drifting (downstream or with advection) propagations. The 
former is a characteristic behavior of RD systems and persists 
while the reaction rate of the activator is dominant (analogues 
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FIG. 5: (color online) Bifurcation diagram showing the branches of 
uniform states (1^0,^0), homoclinic orbits (HO), and traveling waves 
as a function of Da in terms of the maximal value of ^;(^) 
at (a) B = 10.4 and (b) B = 9.6. Solid lines imply linear sta- 
bility, including stability of TW^ to long wave lengths perturba- 
tions 1 9], while Da^ mark the onsets of the linear finite wavenum- 
ber Hopf bifurcation to TW^, respectively. The top inset represents 
the nonuniform states in terms of speed while the large (small) isola 
corresponds to the TW~ family emerging from Da~ dX B — 10.4 
{B — 10.2). The bottom shows HO profiles at locations marked by 
(•); in (a) the two dots mark also the two ends of the HO branch. 
Integration details as in the top panel of Fig.|2]but on distinct periodic 
domains. 

to front dynamics fl3]). While the latter is a consequence 
of low excitation and thus subjected to a nonlinear convec- 
tive instability resulting in a fluid type behavior. Through a 
bifurcation analysis of spatial extended steady states arising 
in a minimal RDA model, we revealed the properties and the 
organization of drifting pulses. Since the results center on 
homoclinic orbits which known to act as organizing centers 
of spatial solutions, qualitative applicability to systems with 
other autocatalytic properties is naturally anticipated. 

Up-to-date only excitable (upstream) solitary waves have 
been observed experimentally in an autocatalytic RDA sys- 
tem [4], nevertheless chemical media operated in cross-flow 
(membrane) tubular reactors |6] or on a rotating disks lIlSll . are 



the most natural setups to confirm our predictions and explore 
technological directions. Moreover, theoretical insights ex- 
plored here can be related to a profound puzzle of large intra- 
cellular particles (organelles) self-organization, in eucaryotic 
cells 1 16]. For example, localized aggregations of myosin-X 
within the filopodia have been observed to propagate bidirec- 
tionally lllTil and from the modeling point of view argued to be 
driven by both diffusion and differential advection [18]. Con- 
sequently, a theoretical framework integrating autocatalytic 
kinetics and distinct transport, is paramount to promoting a 
mechanistic understanding of spatiotemporal trafficking of in- 
tracellular molecular aggregations. 

We thank to N. Gov and M. Naoz for the helpful discussion 
on molecular motors. This work was supported by the US- 
Israel Binational Science Foundation (BSF) and A.Y. was also 
partially supported by the Center for Absorption in Science, 
Israeli Ministry of Immigrant Absorbtion. M.S. is a member 
of the Minerva Center of Nonlinear Dynamics and Complex 
Systems. 



[1] E. Meron, Phys. Rep. 218, 1 (1992); M.C. Cross and RC. Ho- 

henberg, Rev. Mod. Phys. 65, 851 (1993). 
[2] A. Yochehs, E. Knobloch, Y. Xie, Z. Qu, and A. Garfinkel, Eu- 

rophys. Lett. 83, 64005 (2008). 
[3] J.M. Chomaz, Phys. Rev. Lett. 69, 1931 (1992). 
[4] M. K^rn and M. Menzinger, Phys. Rev. E 65, 046202 (2002). 
[5] A. Yochelis and M. Sheintuch, e-print: 

arXiv:nlin.PS/0810.4690; e-print: arXiv:nlin.PS/0902.2688. 
[6] O. Nekhamkina, B.Y Rubinstein, and M. Sheintuch, AIChE J. 

46, 1632 (2000). 

[7] Le (Lewis number), is the ratio of solid- to fluid-phase heat ca- 
pacities, Pe (Peclet number), is the ratio of convective to con- 
ductive enthalpy fluxes, and Da (Damkohler number), is the 
dimensionless rate constant | 6]. 
[8] E. Doedel et ah, AUTO2000: Continuation and bifurcation 
software for ordinary differential equations (with HOMCONT), 
http : //indy. c s . concordia. ca/auto/ 
[9] Temporal stability of TW^ was computed for large periodic 
domains, L = nAc, n > 1, Xc = 27v/kc (until the onset didn't 
change with n), via a standard numerical eigenvalue method us- 
ing Eq. ([T]) in a comoving frame, kc is the critical wavenumber 
at the Hopf onset. 

[10] C. Elphick, E. Meron, J. Rinzel, and E.A. Spiegel, J. Theor. 
Biol. 146, 249 (1990); M. Or-Guil, I.G. Kevrekidis, and M. Bar, 
PhysicaD 135, 154 (2000). 

[11] L.A. Belyakov, Mat. Zametki 28, 910 (1980). 

[12] P Glendinning and C.T. Sparrow, J. Stat. Phys. 35, 645 (1984); 
N.J. Balmforth, Annu. Rev. Fluid Mech. 27, 335 (1995); A.R. 
Champneys et ai, SIAM J. Appl. Dyn. Syst. 6, 663 (2007). 

[13] J. Sneyd, A. LeBeaub, and D. Yule, Physica D 145, 158 (2000). 

[14] V. Yakhnin and M. Menzinger, Chem. Eng. Sci. 57, 4559 
(2002); M. Sheintuch, Y. Smagina, and O. Nekhamkina, Ind. 
Eng. Chem. Res. 41, 2136 (2002). 

[15] Y. Khazan and L.M. Pismen, Phys. Rev. Lett. 75, 4318 (1995). 

[16] M.A. Welte, Curr. Biol. 14, R525 (2004). 

[17] J.S. Berg and R.E. Cheney, Nature Cell. Biol. 4, 246 (2002). 

[18] D.A. Smith and R.M. Simmons, Biophys. J. 80, 45 (2001). 



